Celestin Apprentice 7/Apprentice-Release7.iso/Source Code/Pascal/Applications/NIH Image 1.62b11/src/Dicom.p

Text File | 1996-11-22 | 22KB | 818 lines

{
Dicom.p
  by: Jim Nash, Synergistic Research Systems ([e-mail])
  Reads and decodes the DICOM header so that NIH Image can
  import DICOM images. DICOM (Digital Imaging and Communications
  in Medicine) is a format popular in the medical imaging
  community. This code is in the public domain.
}


unit DICOM;

interface

    uses
        Types, Memory, QuickDraw, Packages, Menus, Events, Fonts, Scrap, ToolUtils, Resources,
        Errors, Palettes, Printing, StandardFile, Folders, TextUtils, Files,
        globals, Utilities, Text, Graphics, Utilities, file2;


    procedure ImportDICOMImages (fname: Str255; RefNum: integer; ImportAll: boolean; {}
                                    function ImportFile (fname: str255; vnum: integer): boolean);

implementation

    const
        dDicomName = 'DICOM dictionary';
        maxElements = 1000;
        elemNameLength = 50;
    type
        DataKind = (kUnknown, kString, kInteger, kLongint, kReal, kUInteger, kULongint);

        DataElement = record
                group, element: integer;
                code: packed array[1..2] of char;
                list: boolean;
                name: string[elemNameLength];
            end;
        ElemArray = array[1..maxElements] of DataElement;
        ElemArrayPtr = ^ElemArray;

        DataDictionary = record
                number: integer;
                elem: ElemArrayPtr;
            end;
        DataDictionaryPtr = ^DataDictionary;
    var
        dictionary: DataDictionaryPtr;
        loaded: boolean;
        mySliceSpacing: real;


{ **************  Utility routines ***************** }

    procedure StringToBase (s: Str255; base: integer; var value: longint);
{converts a string in some base to longint.  Typically}
{base = 2,8,10,16 to represent binary, octal, decimal and hexadecimal}
        var
            ch: char;
            good: boolean;
            len, digit: integer;
            i: longint;
    begin
        i := 1;
        value := 0;
        len := length(s);
        while (i <= len) do begin
            good := true;
            ch := s[i];
            if ch in ['A'..'Z'] then
                digit := ord(ch) - ord('A') + 10
            else if ch in ['a'..'z'] then
                digit := ord(ch) - ord('a') + 10
            else if ch in ['0'..'9'] then
                digit := ord(ch) - ord('0')
            else
                good := false;
            if good then
                value := value * base + digit;
            i := i + 1;
        end;
    end;


    procedure BaseToString (value: longint; base: integer; var s: Str255);
{converts a long integer to a string in any base.  Typically}
{base = 2,8,10,16 to represent binary, octal, decimal and hexadecimal.}
{Ignores the sign bit unless base=10.}
        var
            sign, decimal: boolean;
            digit: integer;
            ch: char;
    begin
        decimal := (base = 10);
        s := '';
        sign := (value < 0);
        if decimal then
            value := abs(value)
        else
            value := BAND(value, $7FFFFFFF);
        if value = 0 then
            s := '0'
        else
            while (value <> 0) do begin
                digit := value mod base;
                value := value div base;
                if (digit >= 0) and (digit <= 9) then
                    ch := chr(digit + ord('0'))
                else
                    ch := chr(digit - 10 + ord('A'));
                s := concat(ch, s);
            end;
        if sign then
            if decimal then
                s := concat('-', s)
            else if s[1] < '2' then
                s[1] := chr(ord(s[1]) + 8)
            else
                s[1] := chr(ord(s[1]) - ord('2') + ord('A'));
    end;


    function htos (i: longint): Str255;
{A convenience function to replace BaseToString (hexadecimal) }
        var
            s: Str255;
    begin
        BaseToString(i, 16, s);
        htos := s;
    end;


    function itos (i: longint): Str255;
{A convenience function to replace NumToString}
        var
            s: Str255;
    begin
        NumToString(i, s);
        itos := s;
    end;


{ **************  DICOM routines ***************** }


    procedure InitDICOM;
        var
            err: integer;
    begin
        dictionary := nil;
        loaded := false;
        DicomInitialized:=true;
    end;



    function FindDicomElement (group, element: integer): integer;
        var
            i, index: integer;
    begin
        index := 0;
        if loaded and (dictionary <> nil) then
            with dictionary^ do
                if (elem <> nil) then begin
                    i := 1;
                    while (group > elem^[i].group) and (i < number) do
                        i := i + 1;
                    if (i <= number) then
                        while (element > elem^[i].element) and (group = elem^[i].group) and (i < number) do
                            i := i + 1;
                    if (i <= number) then
                        if (element = elem^[i].element) and (group = elem^[i].group) then
                            index := i;
                end;
        FindDicomElement := index;
    end;


    procedure InitUserChoice;
{selected elements to list in text window, short form.}
{This a minimum set.  If a user wants more information, they do a full dump.}

        procedure Select (group, element: integer);
            var
                index: integer;
        begin
            with dictionary^ do begin
                index := FindDicomElement(group, element);
                if index > 0 then
                    elem^[index].list := true;
            end;
        end;

    begin
        with dictionary^ do begin
            Select($8, $20);
            Select($8, $30);
            Select($8, $60);
            Select($8, $1030);
            Select($8, $103E);
            Select($8, $1070);

            Select($10, $10);
            Select($10, $20);
            Select($10, $21B0);

            Select($18, $10);
            Select($18, $50);
            Select($18, $88);

            Select($20, $10);
            Select($20, $11);
            Select($20, $12);
            Select($20, $13);

            Select($28, $10);
            Select($28, $11);
            Select($28, $30);
            Select($28, $100);
        end;
    end;


    procedure LoadDataDictionary;
        type
            CharBuf = packed array[0..100000] of char;
            CharBufPtr = ^CharBuf;
        var
            err, refnum, len, i1, i2, n: integer;
            index1, index2, logEOF, count, num, theSize: longint;
            f: text;
            sp: StringPtr;
            str: Str255;
            s1: Str255;
            buf: CharBufPtr;

    begin
        if dictionary = nil then begin
            dictionary := DataDictionaryPtr(NewPtr(sizeof(DataDictionary)));
            if dictionary <> nil then
                dictionary^.elem := ElemArrayPtr(NewPtr(sizeof(ElemArray)));
            loaded := false;
        end;
        if (not loaded) and (dictionary <> nil) then
            with dictionary^ do begin
                err := HSetVol(nil, StartupSpec.vRefNum, StartupSpec.parID);
                err := FSOpen(dDicomName, 0, refnum);     {check that file is present}
                if (err = 0) and (elem <> nil) then begin
                    err := GetEOF(refnum, logEOF);
                    buf := CharBufPtr(NewPtr(logEOF + 10));
                    if (buf <> nil) then begin
                        loaded := true;
                        number := 0;
                        count := logEOF;
                        err := FSRead(refnum, count, ptr(buf));
                        err := FSClose(refnum);
                        index1 := 0;
                        repeat
                            index2 := index1;
                            str := '';
                            while (buf^[index2] <> cr) and (index2 < logEOF) and (length(str) < 255) do begin
                                str := concat(str, buf^[index2]);
                                index2 := index2 + 1;
                            end;
                            index1 := index2 + 1;
                            len := length(str);
                            if len > 0 then
                                if str[1] = '{' then begin
                                    number := number + 1;
                                    if (number mod 10) = 0 then
                                        ShowAnimatedWatch;
                                    with elem^[number] do begin
                                        list := false;

                                        i1 := pos('x', str);
                                        s1 := copy(str, i1 + 1, 4);
                                        StringToBase(s1, 16, num);
                                        group := num;
                                        str := copy(str, i1 + 6, length(str)-(i1 + 6));

                                        i1 := pos('x', str);
                                        s1 := copy(str, i1 + 1, 4);
                                        StringToBase(s1, 16, num);
                                        element := num;
                                        str := copy(str, i1 + 6, length(str)-(i1 + 6));

                                        i1 := pos('''', str);
                                        if length(str) >= (i1 + 2) then begin
                                            code[1] := str[i1 + 1];
                                            code[2] := str[i1 + 2];
                                            str := copy(str, i1 + 5, length(str)-(i1 + 5));
                                        end
                                        else
                                            str := '';

                                        i1 := pos('"', str);
                                        if i1 > 0 then
                                            str[i1] := ' ';
                                        i2 := pos('"', str);
                                        if i2 = 0 then
                                            number := number - 1
                                        else begin
                                            n := i2 - i1 - 1;
                                            if n > elemNameLength then
                                                n := elemNameLength;
                                            name := copy(str, i1 + 1, n);
                                        end;
                                    end;
                                end;
                        until (index1 >= logEOF);
                    end;
                    DisposePtr(ptr(buf));
                end;
                InitUserChoice;
            end;
    end;
{$R+}


    function GetDataKind (index: integer): DataKind;
        var
            kind: DataKind;
    begin
        kind := kUnknown;
        if (dictionary <> nil) and (index > 0) then
            with dictionary^.elem^[index] do begin
                if (code = 'AE') or (code = 'AS') or (code = 'CS') or (code = 'DA') or (code = 'DS') then
                    kind := kString
                else if (code = 'DT') or (code = 'IS') or (code = 'LO') or (code = 'LT') or (code = 'PN') then
                    kind := kString
                else if (code = 'SH') or (code = 'ST') or (code = 'TM') or (code = 'UI') then
                    kind := kString
                else if (code = 'SS') then
                    kind := kInteger
                else if (code = 'SL') then
                    kind := kLongint
                else if (code = 'US') then
                    kind := kUInteger
                else if (code = 'UL') then
                    kind := kULongint;
            end;
        GetDataKind := kind;
    end;

    procedure ImportDICOMImages (fname: Str255; RefNum: integer; ImportAll: boolean; {}
                                    function ImportFile (fname: str255; vnum: integer): boolean);
        var
            enable_text, enable_open_text, first_image, sw, listAll, UseFixedScale: boolean;
            ImageNumber:integer;
            myIntercept, myScale: extended;

        function GetDICOMParams (fname: Str255; vNum: integer): integer;
            const
                id_offset = 128;            {location of "DICM"}
                firstDicomElement = 132;    {first element}
                maxbuf = 20000;
            type
                name4 = packed array[1..4] of char;
                name4ptr = ^name4;
            var
                open_sw, done, window_sw: boolean;
                f, err, index, len: integer;
                groupWord, elementWord, lastGroup, FirstElement: integer;
                height, width, bits_alloc, bits_stored, high_bit, representation, offset, bitsAllocated: integer;
                seriesMin, seriesMax: integer;           {intensity range}
                scale, aspect, units: Str255;       {spatial}
                s, imgNumString, sliceSpacingStg, rescaleInterceptStg, rescaleSlopeStg: Str255;
                buflen, elementLength, groupLength: longint;
                buf: packed array[0..maxbuf] of byte;
                kind: DataKind;
                dictionaryIndex: integer;
                xStr,yStr, vr:str255;

            procedure MyWriteElement (str: Str255);
                const
                    spaces = '                                                                                               ';
                    padWidth = 4;
                    nameWidth = 32;
                var
                    s1, s2: str255;

                function pad (s: Str255): Str255;
                    const
                        width = 4;
                begin
                    while length(s) < width do
                        s := concat(' ', s);
                    pad := s;
                end;

            begin
                with dictionary^.elem^[dictionaryIndex] do begin
                    s2 := name;
                    if listAll then begin
                        s1 := concat('(', pad(htos(groupWord)), ',', pad(htos(elementWord)), ')  (', pad(itos(elementLength)), ')');
                        s2 := copy(concat(name, spaces), 1, nameWidth);
                        s2 := concat(s1, '  ', code[1], code[2], '  ', s2);
                    end;
                    str := concat(s2, ':  ', str);
                    if enable_text then begin
                        if groupWord <> lastGroup then
                            InsertText('', true);
                        lastGroup := groupWord;
                        InsertText(str, true);
                    end;
                end;
            end;

            function GetInteger (index: integer): integer;
                var
                    i: integer;
            begin
                i := buf[index] + $100 * buf[index + 1];
                GetInteger := i;
            end;

            function GetLongint (index: integer): longint;
                var
                    i: longint;
            begin
                i := Ord4(buf[index]) + $100 * (buf[index + 1] + $100 * (buf[index + 2] + $100 * buf[index + 3]));
                GetLongint := i;
            end;

            function GetLength (i: integer): longint;
            begin
                vr[1] := chr(buf[i]);
                vr[2] := chr(buf[i + 1]);
                if vr[2] < 'A' then {implecit vr with 32-bit length}
                    GetLength := Ord4(buf[i]) + $100 * (buf[i+1] + $100 * (buf[i+2] + $100 * buf[i+3]))
                else if (vr = 'OB') or (vr = 'OW') or (vr = 'SQ') then begin  {explicit VR with 32-bit length}
                    i := i + 4;  {skip 2 byte string and 2 reserved bytes}
                    index := index + 4;
                    GetLength := Ord4(buf[i]) + $100 * (buf[i+1] + $100 * (buf[i+2] + $100 * buf[i+3]))
                end else  {explicit VR with 16-bit length}
                    GetLength := Ord4(buf[i+2]) + $100 * (buf[i+3]);
            end;

            function GetUInteger (index: integer): longint;
                var
                    i: integer;
            begin
                i := buf[index] + $100 * buf[index + 1];
                GetUInteger := BAND(i, $FFFF);
            end;

            function GetULongint (index: integer): longint;
            {does not correctly report numbers > $7FFFFFFF}
            begin
                GetULongint := GetLongint(index);
            end;

            function GetString (index: integer): Str255;
                var
                    i: integer;
                    s: Str255;
            begin
                s := '';
                for i := 0 to elementLength - 1 do
                    s := concat(s, chr(buf[index + i]));
                GetString := s;
            end;


            function htos2(i: LongInt): str255;
            {Converts an integer to hex using a fixed field width of 6}
            var
                s: str255;
            begin
                s := htos(i);
                while length(s) < 6 do
                    s := concat(' ', s);
                htos2 := s;
            end;


            procedure DoByteSwap (var i: LongInt);
                var
                    a: ostype;
                    c: char;
            begin
                a := ostype(i);
                c := a[1];
                a[1] := a[2];
                a[2] := c;
                c := a[3];
                a[3] := a[4];
                a[4] := c;
                i := LongInt(a)
            end;


    procedure Swap(var i: LongInt);
        var
            a: ostype;
            c: char;
    begin
        a := ostype(i);
        a[1] := a[3];
        a[2] := a[4];
        a[3] := chr(0);
        a[4] := chr(0);
        i := LongInt(a)
    end;


            procedure GetNextElement;
                var
                    i: longint;
                    str: Str255;
            begin
                if index = 0 then
                    index := firstElement
                else
                    index := index + elementLength;
                if (index < 0) or (index >= buflen) then
                    exit(GetNextElement);
                groupWord := GetInteger(index);
                elementWord := GetInteger(index + 2);
                elementLength := GetLength(index + 4);
                if elementLength = 13 then {hack needed to read some GE files}
                    elementLength := 10;
                if ControlKeyDown then
                    InsertText(stringOf(index:6, htos2(groupWord), htos2(elementWord), htos2(elementLength)), true);
                index := index + 8;
                dictionaryIndex := FindDicomElement(groupWord, elementWord);
                if dictionaryIndex > 0 then
                    with dictionary^ do
                        if elem^[dictionaryIndex].list or listAll then
                            with elem^[dictionaryIndex] do begin
                                kind := GetDataKind(dictionaryIndex);
                                case kind of
                                    kString:
                                        str := GetString(index);
                                    kInteger:
                                        str := itos(GetInteger(index));
                                    kLongint:
                                        str := itos(GetLongint(index));
                                    kUInteger:
                                        str := itos(GetLongint(index));
                                    kULongint:
                                        str := itos(GetULongint(index));
                                    otherwise
                                        str := 'unknown format';
                                end;
                                MyWriteElement(str);
                            end;
            end;

            function IsElement (group, element: integer): boolean;
            begin
                IsElement := (group = groupWord) and (element = elementWord);
            end;

            procedure TestError (err1: integer; str: Str255);
                var
                    str1: Str255;
            begin
                err := err1;
                if err1 <> 0 then begin
                    if err <> 1 then
                        str := concat(str, ' - error ', itos(err));
                    PutMessage(str);
                    if open_sw then
                        err1 := fsclose(f);
                    GetDICOMParams := err;
                    exit(GetDICOMParams);
                end;
            end;

            procedure OpenDicomTextWindow;
                var
                    width, height: integer;
            begin
                if listAll then begin
                    width := 500;
                    height := 400;
                end
                else begin
                    width := 350;
                    height := 300;
                end;
                if enable_open_text then
                    window_sw := MakeNewTextWindow(concat(fname, ' header'), width, height);
                CurrentFontID := monaco;
                CurrentSize := 9;
                ChangeFontOrSize;
                enable_open_text := false;
                if enable_text then begin
                    if loaded then
                        InsertText('Selected fields from the DICOM file header', true)
                    else begin
                        InsertText(concat('Can''t find file: ', dDicomName, '.'), true);
                        InsertText('', true);
                        InsertText('This file is required to decode the DICOM header. It is available from: ftp://zippy.nimh.nih.gov/pub/nih-image/documents/dicom-dict.hqx. It must be located in the same folder as NIH Image or in the System folder.', true)
                    end;
                    InsertText('', true);
                end;
            end;

        begin
            vr :='--';
            err := 0;
            buflen := maxbuf + 1;
            open_sw := false;
            TestError(fsopen(fname, vNum, f), 'Open');
            open_sw := true;
            TestError(FSRead(f, buflen, @buf), 'Read');
            TestError(fsclose(f), 'Close');
            if name4ptr(longint(@buf) + id_offset)^ = 'DICM' then
                FirstElement:=FirstDicomElement
            else if name4ptr(longint(@buf))^ = 'DICM' then
                FirstElement:=4
            else
                FirstElement:=0; {TestError(1, 'This is not a DICOM file.');}
            OpenDicomTextWindow;
            sliceSpacingStg := '1.0';
            rescaleInterceptStg := '0.0';
            rescaleSlopeStg := '1.0';
            imgNumString := '';
            height := -1;
            width := -1;
            offset := -1;
            index := 0;
            lastGroup := $8;
            done := false;
            scale := '0.0';
            aspect := '1.0';
            representation := 0;
            bitsAllocated := 16;
            seriesMin := 0;
            seriesMax := 0;
            units := '';
            repeat
                GetNextElement;
                if (index < 0) or (index >= buflen) then
                    leave;
                if (elementWord = 0) and (elementLength = 0) then
                    leave;

                if IsElement($18, $88) then
                    sliceSpacingStg := GetString(index)
                else if IsElement($20, $13) then
                    imgNumString := GetString(index)
                else if IsElement($28, $10) then
                    height := GetInteger(index)
                else if IsElement($28, $11) then
                    width := GetInteger(index)
                else if IsElement($28, $30) then
                    scale := GetString(index)
                else if IsElement($28, $34) then
                    aspect := GetString(index)
                else if IsElement($28, $100) then
                    bitsAllocated := GetInteger(index)
                else if IsElement($28, $103) then
                    representation := GetInteger(index)
                else if IsElement($28, $108) then
                    seriesMin := GetInteger(index)
                else if IsElement($28, $109) then
                    seriesMax := GetInteger(index)
                else if IsElement($28, $1052) then
                    rescaleInterceptStg := GetString(index)
                else if IsElement($28, $1053) then
                    rescaleSlopeStg := GetString(index)
                else if IsElement($7FE0, $10) then begin
                    offset := index;
                    done := true;
                end;
                if CommandPeriod then
                    listAll := false;
            until done;
            if (width = -1) or (height = -1) or (offset = -1) then begin
                InsertText('', true);
                InsertText('Unable to decode DICOM header.', true);
                GetDICOMParams := -1;
                exit(GetDICOMParams)
            end;

            {Image dimension information}
            ImportCustomWidth := width;
            ImportCustomHeight := height;
            ImportCustomOffset := offset;
            ImportSwapBytes := true; {(representation = 1);}

            {Intensity information}
            if bitsAllocated = 8 then begin
                ImportCustomDepth := EightBits;
                if ImageNumber=1 then
                    InsertText('', true);
                GetDICOMParams := err;
                exit(GetDICOMParams);
            end else
                ImportCustomDepth := SixteenBitsSigned;
            if not ((seriesMin = 0) and (seriesMax = 0)) then begin
                ImportAutoScale:=false;
                ImportMin:=seriesMin;
                ImportMax:=seriesMax;
            end else if (ImageNumber>1) and UseFixedScale then begin
                ImportAutoScale:=false;
                ImportMin:=info^.CurrentMin;
                ImportMax:=info^.CurrentMax;
            end else
                ImportAutoScale:=true;;

            {convert from scaled units to independent units}
            myIntercept := StringToReal(rescaleInterceptStg);
            myScale := StringToReal(rescaleSlopeStg);

            {Spatial scale information}
            with Info^ do begin
                PixelAspectRatio := StringToReal(aspect);
                xScale := 1;
                yScale := 1;
                zScale := 1.0 / StringToReal(sliceSpacingStg);

                xUnit := '';
                SpatiallyCalibrated := false;
                if scale <> '' then begin
                    xStr:=copy(scale, pos('\', scale) + 1, length(scale) - pos('\', scale));
                    xScale := StringToReal(xStr);
                    yStr:=copy(scale, 1, pos('\', scale) - 1);
                    yScale := StringToReal(yStr);
                    xUnit := 'mm';
                    SpatiallyCalibrated := (xScale <> 0.0) and (yScale <> 0.0);
                    if SpatiallyCalibrated then begin
                        xScale := 1.0 / xScale;
                        yScale := 1.0 / yScale;
                    end;
                end;
            end; {with}
            if ImageNumber=1 then
                InsertText('', true);
            GetDICOMParams := err;
        end;


        procedure UpdateCoefficients;
        {Scale coefficients given Dicom Rescale Intercept and Rescale Slope}
        begin
            with info^ do begin
                info^.Coefficient[1] := myIntercept + myScale * info^.Coefficient[1];
                info^.Coefficient[2] := myScale * info^.Coefficient[2];
                fit := StraightLine;
                GenerateValues;
            end;
        end;


        procedure ImportAllDicomFiles (RefNum: integer);
            var
                OpenedOK: boolean;
                index: integer;
                name: Str255;
                ftype: OSType;
                err: OSErr;
                PB: HParamBlockRec;
        begin
            index := 0;
            while true do begin
                index := index + 1;
                with PB do begin
                    ioCompletion := nil;
                    ioNamePtr := @name;
                    ioVRefNum := RefNum;
                    ioVersNum := 0;
                    ioFDirIndex := index;
                    err := PBGetFInfoSync(@PB);
                    if err = fnfErr then
                        exit(ImportAllDicomFiles);
                    ftype := ioFlFndrInfo.fdType;
                end;

                if GetDICOMParams(name, RefNum) <> 0 then
                    exit(ImportAllDicomFiles);
                WhatToImport := ImportCustom;
                if not ImportFile(name, RefNum) then
                    exit(ImportAllDicomFiles);
                WhatToImport := ImportDicom;
                if (myIntercept <> 0.0) or (myScale <> 1.0) then
                    UpdateCoefficients;
                with info^ do InsertText(StringOf(ImageNumber:3, ': "', title, '", min=', CurrentMin:1, ', max=', CurrentMax:1), true);
                ImageNumber:=ImageNumber+1;
                enable_text := false;        {text saved for first image only}
                first_image := false;
                if CommandPeriod then begin
                    beep;
                    exit(ImportAllDicomFiles);
                end;
            end;
        end;

    begin        {ImportDICOMImages}
        if not DicomInitialized then
            InitDICOM;
        listAll := OptionKeyDown or OptionKeyWasDown;
        enable_open_text := true;
        enable_text := true;
        first_image := true;
        ImageNumber:=1;
        UseFixedScale:=ShiftKeyDown;
        LoadDataDictionary;
        ImportingDicom := true;
        if ImportAll then
            ImportAllDicomFiles(RefNum)
        else begin
            if GetDICOMParams(fname, RefNum) <> 0 then
                exit(ImportDICOMImages);
            WhatToImport := ImportCustom;
            if ImportFile(fname, RefNum) then
                if (myIntercept <> 0.0) or (myScale <> 1.0) then
                    UpdateCoefficients;
            WhatToImport := ImportDicom;
            with info^ do InsertText(StringOf('file "', title, '", min=', CurrentMin:1, ', max=', CurrentMax:1), true);
        end;
        ImportingDicom := false;
    end;

end.